home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / RANDPORT.C < prev    next >
C/C++ Source or Header  |  1995-09-07  |  6KB  |  238 lines
  1. /* Listing 2
  2.     rand_por[t].c
  3.     a portable and reasonably fast random number generator
  4.         multiplicative random number generator
  5.     see
  6.        L'Ecuyer - Comm. of the ACM, Oct. 1990, vol. 33.
  7.        Numerical Recipes in C, 2nd edition, pp.  278-86
  8.        L'Ecuyer and Cote, ACM Transactions on Mathematical Software,
  9.          March 1991
  10.        Russian peasant algorithm -- Knuth, vol. II, pp.  442-43
  11.     Copyright (c) 1994, 1995 by Gerald P. Dwyer, Jr.     */
  12.  
  13. /* ------------------- */
  14. /* FUNCTION PROTOTYPES */
  15. /* ------------------- */
  16. # undef F
  17. # if defined(__STDC__) || defined(__PROTO__)
  18. #    define  F( P )  P
  19. # else
  20. #    define  F( P )  ()
  21. # endif
  22.  
  23. /* INDENT OFF */
  24. extern    long    genr_rand_port F((long));
  25. extern    long    get_init_rand_port F((void));
  26. extern    long    init_rand_port F((long));
  27. extern    long    mult_mod F((long, long, long));
  28. extern    long    rand_port F((void));
  29. extern    double    rand_rect_port F((void));
  30. extern    long    skip_ahead F((long, long, long, long));
  31.  
  32. # undef F
  33. /* INDENT ON */
  34.  
  35. #include    <time.h>
  36. #include    <stdlib.h>
  37. #include    <limits.h>
  38. #include    <assert.h>
  39.  
  40. #define    TESTING
  41. #undef    TESTING
  42.  
  43. #define  MOD    2147483647L     /* modulus for generator */
  44. #define  MULT        41358L     /* multiplier            */
  45.                           /* modulus = mult*quotient + remainder */
  46. #define  Q           51924L     /* int(modulus / multiplier) */
  47. #define  R           10855L     /* remainder                 */
  48. #define  MAX_VALUE    (MOD-1)
  49.  
  50. #define  EXP_VAL 1285562981L    /* value for 10,000th draw */
  51.  
  52. #define  IMPOSSIBLE_RAND (-1)
  53. #define  STARTUP_RANDS   16     /* throw away this number of
  54.                           initial random numbers */
  55.  
  56. static long rand_num = IMPOSSIBLE_RAND ;
  57.  
  58. /* initialize random number generator with seed */
  59. long init_rand_port(long seed)
  60. {
  61.     extern long rand_num ;
  62.     int i ;
  63.  
  64.     if (seed < 1 || seed > MAX_VALUE)        /* if seed out of range */
  65.         seed = get_init_rand_port() ;            /* get seed */
  66.  
  67.     rand_num = seed ;
  68.     for (i = 0; i < STARTUP_RANDS; i++)    /* and throw away */
  69.         rand_num = genr_rand_port(rand_num) ;    /* some initial ones */
  70.  
  71.     return seed ;
  72. }
  73.  
  74.  
  75. /* get a long initial seed for generator
  76.     assumes that rand returns a short integer */
  77. long get_init_rand_port(void)
  78. {
  79.     long seed ;
  80.  
  81.     srand((unsigned int)time(NULL)) ;    /* initialize system generator */
  82.     do {
  83.         seed  = ((long)rand())*rand() ;
  84.         seed += ((long)rand())*rand() ;
  85.     } while (seed > MAX_VALUE) ;
  86.  
  87.     assert (seed > 0) ;
  88.  
  89.     return seed ;
  90. }
  91.  
  92.  
  93. /* generate the next value in sequence from generator
  94.     uses approximate factoring
  95.     residue = (a * x) mod modulus
  96.             = a*x - [(a*x)/modulus]*modulus
  97.     where
  98.         [(a*x)/modulus] = integer less than or equal to (a*x)/modulus
  99.     approximate factoring avoids overflow associated with a*x and
  100.         uses equivalence of above with
  101.     residue = a * (x - q * k) - r * k + (k-k1) * modulus
  102.     where
  103.         modulus = a * q + r
  104.         q  = [modulus/a]
  105.         k  = [x/q]        (= [ax/aq])
  106.         k1 = [a*x/modulus]
  107.     assumes
  108.         a, m > 0
  109.         0 < init_rand < modulus
  110.         a * a <= modulus
  111.         [a*x/a*q]-[a*x/modulus] <= 1
  112.             (for only one addition of modulus below) */
  113. long genr_rand_port(long init_rand)
  114. {
  115.     long k, residue ;
  116.  
  117.     k = init_rand / Q ;
  118.     residue = MULT * (init_rand - Q * k) - R * k ;
  119.     if (residue < 0)
  120.         residue += MOD ;
  121.  
  122.     assert(residue >= 1 && residue <= MAX_VALUE) ;
  123.  
  124.     return residue ;
  125. }
  126.  
  127.  
  128. /* get a random number */
  129. long rand_port(void)
  130. {
  131.     extern long rand_num ;
  132.  
  133.                 /* if not initialized, do it now */
  134.     if (rand_num == IMPOSSIBLE_RAND) {
  135.         rand_num = 1 ;
  136.         init_rand_port(rand_num) ;
  137.     }
  138.  
  139.     rand_num = genr_rand_port(rand_num) ;
  140.  
  141.     return rand_num ;
  142. }
  143.  
  144.  
  145. /* generates a value on (0,1) with mean of .5
  146.     range of values is [1/(MAX_VALUE+1), MAX_VALUE/(MAX_VALUE+1)]
  147.     to get [0,1], use (double)(rand_port()-1)/(double)(MAX_VALUE-1) */
  148. double rand_rect_port(void)
  149. {
  150.     return (double)rand_port()/(double)(MAX_VALUE+1) ;
  151. }
  152.  
  153.  
  154. /* skip ahead in recursion
  155.     residue = (a^skip * init) mod modulus
  156.     Use Russian peasant algorithm          */
  157. long skip_ahead(long a, long init_rand, long modulus, long skip)
  158. {
  159.     long residue = 1 ;
  160.  
  161.     if (init_rand < 1 || init_rand > modulus-1 || skip < 0)
  162.         return -1 ;
  163.  
  164.     while (skip > 0) {
  165.         if (skip % 2)
  166.             residue = mult_mod(a, residue, modulus) ;
  167.         a = mult_mod(a, a, modulus) ;
  168.         skip >>= 1 ;
  169.     }
  170.     residue = mult_mod(residue, init_rand, modulus) ;
  171.  
  172.     assert(residue >= 1 && residue <= modulus-1) ;
  173.  
  174.     return residue ;
  175. }
  176.  
  177.  
  178. /* calculate residue = (a * x) mod modulus for arbitrary a and x
  179.         without overflow
  180.     assume 0 < a < modulus and 0 < x < modulus
  181.     use Russian peasant algorithm followed by approximate factoring */
  182. long mult_mod(long a, long x, long modulus)
  183. {
  184.     long q, r, k, residue ;
  185.  
  186.     residue = -modulus ;        /* to avoid overflow on addition */
  187.  
  188.     while (a > SHRT_MAX) {    /* use Russian Peasant to reduce a */
  189.         if (a % 2) {
  190.             residue += x ;
  191.             if (residue > 0)
  192.                 residue -= modulus ;
  193.         }
  194.         x += (x - modulus) ;
  195.         if (x < 0)
  196.             x += modulus ;
  197.         a >>= 1 ;
  198.     }
  199.                     /* now apply approximate factoring to a
  200.                         and compute (a * x) mod modulus */
  201.     q = modulus / a ;
  202.     r = modulus - a * q ;
  203.     k = x / q ;
  204.     x = a * (x - q * k) - r * k ;
  205.     while (x < 0)
  206.         x += modulus ;
  207.                     /* add result to residue and take mod */
  208.     residue += x ;
  209.     if (residue < 0)    /* undo initial subtraction if necessary */
  210.         residue += modulus ;
  211.  
  212.     assert(residue >= 1 && residue <= modulus-1) ;
  213.  
  214.     return residue ;
  215. }
  216.  
  217.  
  218. # if defined(TESTING)
  219. /* Test the generator */
  220. #include    <stdio.h>
  221. void main(void)
  222. {
  223.     long seed ;
  224.     int i ;
  225.  
  226.     seed = init_rand_port(1) ;
  227.     printf("Seed for random number generator is %ld\n", seed) ;
  228.     i = STARTUP_RANDS ;   /* threw away STARTUP_RANDS */
  229.     do {
  230.         rand_port() ;
  231.         i++ ;
  232.     } while (i < 9999) ;
  233.  
  234.     printf("On draw 10000, random number should be %ld\n", EXP_VAL) ;
  235.     printf("On draw %d, random number is        %ld\n", i+1, rand_port()) ;
  236. }
  237. # endif    /* TESTING */
  238.